Effects of intrinsic fluctuations in a prototypical chemical oscillator: metastability and switching 
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Intrinsic or demographic noise has been shown to play an important role in the dynamics of a variety of sys- 
tems including predator-prey populations, intracellular biochemical reactions, and oscillatory chemical reaction 
systems, and is known to give rise to oscillations and pattern formation well outside the parameter range pre- 
dicted by standard mean-field analysis. Initially motivated by an experimental model of cells and tissues where 
the cells are represented by chemical reagents isolated in emulsion droplets, we study the stochastic Brusse- 
lator, a simple activator-inhibitor chemical reaction model. Our work builds on the results of recent studies 
and looks to understand the role played by intrinsic fluctuations when the timescale of the inhibitor species is 
fast compared to that of the activator. In this limit, we observe a noise induced switching between small and 
large amplitude oscillations that persists for large system sizes (N), and deep into the non oscillatory part of the 
mean-field phase diagram. We obtain a scaling relation for the first passage times between the two oscillating 
states . From our scaling function, we show that the first passage times have a well defined form in the large A' 
limit. Thus in the limit of small noise and large timescale separation a careful treatment of the noise will lead to 
a set of non- trivial deterministic equations different from those obtained from the standard mean-field limit. 

PACS numbers: 05.40.-a, 02.50.Ey, 82.40.Bj 



I. INTRODUCTION 

In recent years systems with demographic or intrinsic fluc- 
tuations have gained an increasing amount of attention. The 
fluctuations in these systems arise from the stochastic nature 
of discrete reactions between a finite number of elements. 
Careful treatment of the intrinsic noise has lead to the discov- 
ery of many interesting behaviors not found in the traditional 
mean-field treatment of such problems. McKane and Newman 
showed that demographic noise can induce noisy oscillations, 
or quasicycles, where mean-field predicts a stable stationary 
state in a simple predator-prey model iH. Butler and Golden- 
feld then extended this to show that a similar mechanism can 
lead to noisy Turing patterns, quasipattems, well outside the 
region of parameter space predicted by standard Turing anal- 
ysis Butler's work gives a possible solution to the fine 
tuning problem for pattern formation and provides a possible 
explanation for the ubiquity of Turing-like patterns observed 
in nature. The eff'ects of intrinsic stochasticity has also been 
looked at in models of chemical oscillators and cell signal- 
ing and intracellular biochemical reactions, all of which show 
interesting non mean-field behavior llH-l^ . 

In this paper we are interested in systems in which the dy- 
namics of one variable is fast with respect to the other vari- 
ables in the system. This timescale separation can lead to ex- 
citability, meaning that if the system is sufficiently perturbed 
from a stable fixed point, it will go on a large excursion in 
phase space before returning to the fixed point. Examples of 
excitable systems include neuronal systems, and chemical os- 
cillators such as the Belouzov-Zhabotinsky reaction ||7|. In 
recent years, it has been recognized that, in the presence of 
a large time scale separation, even small noise can lead to 
dynamical patterns that are absent in the mean-field model 
Isl- UOtl . It is known that in systems with well-separated time 
scales, a singular Hopf bifurcation can occur leading to a large 
jump in oscillation amplitude and a transition to relaxation os- 
cillations 1 1 111 . The appearance of stochastic resonance has 



been demonstrated in a Brusselator close to a Hopf bifurca- 
tion lfl2ll . This stochastic resonance occurs when the dynam- 
ics of the activator species is fast compared to the inhibitor 
In this paper, we study the efifects of intrinsic fluctuations in 
the opposite limit where the inhibitor species has the fast dy- 
namics. Our results show that , in this regime, the Brusselator 
switches between small and large amplitude oscillations with 
switching times that remain finite in the limit of large systems 
(small intrinsic noise). The switchability is a property of the 
intrinsically noisy system and is not present in the meanfield 
model. 

The remainder of this paper is structured as follows. In 
the next section we will present our version of the Brussela- 
tor model starting with the discrete reactions. Then we will 
briefly describe the behavior of the Brusselator in the mean- 
field limit. This section will conclude with a short discussion 
of the van Kampen system size expansion of the Brussela- 
tor model, which was included for completeness. For a more 
detailed discussion one should see references (jH [HI [l3l- In 
section Hn] we will present the results of our simulations and 
analysis, all of which will be summarized in section HV] 



II. THE BRUSSELATOR MODEL 



The Brusselator, put forth by Ilya Prigogine in the 1960s, 
is a prototypical model of an autocatalytic chemical reaction 
system showing limit cycle behavior iflSll . This model has 
been studied extensively over the past several decades includ- 
ing works on the mean-field limit, subject to external fluctu- 
ations, and more recently with noise due to intrinsic or de- 
mographic noise iH IH [l2 The Brusselator model 
consists of two species, an activator X and inhibitor Y, which 



2 




2 4 6 8 10 



b 

FIG. 1 : Phase diagram of the mean-field Bmsselator. The Brussela- 
tor has a stable Fixed point in regions / and // and a stable limit cycle 
about an unstable fixed point in regions /// and IV. The eigenvalues 
of the Jacobian matrix are complex in regions // and ///, making 
the fixed point a stable, or unstable, spiral, while the eigenvalues in 
regions / and IV are real. 

undergo the set of four reactions 



2X + Y 3X. (1) 

We assume the system is coupled to some external bath where 
X particles are fed into the system with rate and leave with 
rate 1, which set the size and timescale for our system respec- 
tively, and thus the system is maintained out of equilibrium. 
The parameter b is the rate of exchange from X to Y, while c 
is the rate of the autocatalytic reaction converting Y back to 
X. As will be shown below, c also determines the separation 
of timescales between the dynamics of the two species X and 
Y. 



A. Mean-Field - Linear Stability Analysis 

Using the law of mass action, we can write down a set of 
mean-field rate equations for the number density of the two 
species 

i = 1 -(1 +b)x + cx^y (2) 
y — bx — cx^y (3) 

where x = X/N and y = Y/N. These equations have a unique 
fixed point at x* = l,y* - b/c that is linearly stable for b < 
c + I and becomes unstable, giving rise to a stable limit cycle 
as the control parameter, b, is varied past the critical value 
br = c +1. 



The phase diagram, to linear order, is shown in Fig. [T] In 
regions / and // the fixed point is stable, and in /// and IV 
the fixed point is unstable and surrounded by a stable limit 
cycle. The eigenvalues of the Jacobian matrix in regions // 
and /// are complex making the fixed point an unstable (or 
stable) spiral. 



B. Mean-Field - Beyond Linear Stability 

In the presence of a large time scale separation many os- 
cillatory systems will become excitable, meaning that small 
perturbations from the stable fixed point or small amplitude 
limit cycle can send the system on a large amplitude excur- 
sion through the phase space. This time scale separation can 
be written explicitly by the change of variables x' - x and 
y' - cy in Eqs.|2]and|3] Carrying out this substitution gives 

X' = 1 -(1 +b)x' +x'y (4) 
Tyy' - bx' - x'^y' (5) 

where t,. = c"' is the time scale of the variable y. In the large c 
limit, the inhibitor has the fast dynamics whereas as c — » 0, the 
activator becomes the fast variable. Within mean field theory, 
there is a marked change in behavior as the ratio of timescales 
is varied, as illustrated in Fig. |2] The vector plots of the 
phase space along with the nullclines, at c = 0.01, c =1.0 and 
c - 9.0 and b - 0.9, b - 1.9, and b - 9.9 below the bifurca- 
tion threshold, show that as c — > 0, the fixed point approaches 
the peak of the X-variable nullcline, and becomes pinned at 
that point. As c — > c», on the other hand, the segments of the 
X and Y nullclines that lie to the right of the peak, approach 
each other and merge converting the single fixed point to a line 
of fixed points. Sample trajectories shown in Fig. |2]demon- 
strate that when the time scales are comparable, the trajectory 
will simply spiral into the fixed point. When c is large, the 
trajectory will escape to the nullcline of the slow variable (X), 
which it will follow until it quickly jumps to the other slow 
branch following that and ultimately spiraling into the fixed 
point along the incipient small limit cycle. The second slow 
branch is almost superposed on the Y nullcline, and becomes 
the line of fixed points in the c ^ oo limit. When c is small, 
as in Fig. |2^, the trajectory escapes along aX - -7 path, then 
quickly jumps to the slow branch following it to the stable 
fixed point. 

In all three cases, as the control parameter b is tuned past 
the bifurcation point b^ = c+l the fixed point becomes unsta- 
ble and gives rise to a stable limit cycle with an amplitude that 
grows as y/b - be near be. In the regime where c » 1 the am- 
plitude of the limit cycle will jump sharply as it is increased to 
the point where the trajectories can escape. At this point, the 
large excursion becomes the stable limit cycle. Fig. |3]shows 
the amplitude of oscillation as a function of the distance from 
the bifurcations, 6 = b-be. It should be emphasized that in the 
mean field description the Brusselator is alway a monostable 
system being either at a fixed point, a small amplitude limit 
cycle or a large amplitude limit cycle. As we show below, the 
presence of noise induces metastability and switching. 
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FIG. 2: Phase space vector plot of Brusselator equations with nullclines(dashed lines) and example trajectory(solid line) starting near the fixed 
point in three different timescale regiemes; c << l(left), c = l.O(center), and c » l.O(right). 




FIG. 3: Amplitude of oscillation as b is varied past the bifurcation 
point when x and y have the same time scale, c = I.O (□), and when 
y is the fast variable, c = 9.0 (O)- Amplitudes were taken from 
numerical solutions to Eqs.|2]and[3] 



C. Master Equation Formulation 

The rate equations analyzed in the previous section are of 
a mean field nature |(3|, which ignores the demographic noise 
present in the intrinsically stochastic process described by Eq. 
[1] In order to account for the stochasticity, the "random walk" 
in chemical space is written down in the form of a Master 
equation IlKlsKl^l. The relevant equation for the Brusselator 



+(6;-l)r2(n,,nv) 

+(e;e;-l)r3(n„«,) 

+(e;< - l)r4(n,, UyWin,, n,, f). (6) 



Here Pin.^, riy, f) is the time dependent probability of finding 
the system in the state (n v, «,,) at time t and is a step operator 
that acts on functions of tiy^in^) by raising or lowing n^iiiy) by 



one, as follows: 

etf{nx,ny) = f{n^+l,ny) 
£yf{nx,ny) = f{n^,ny± 1). 

and the transition rates, r,(n i, «,,) are given by 



(7) 
(8) 



Ty{n^,ny) = 

T2{nx,ny) = 

Ti{n^,ny) = 

r4(«^,nv) = 



n.v 

w ■ 



(9) 



While the full Master equation cannot be exactly solved, 
significant insight can be gained by carrying out a system size 
expansion. We first assume that the number density of each 
species is given by the average density plus fluctuations of 
order A^^'^^ 



-± ^ x{t) + N-''^Ut) 



(10) 
(11) 



It can be shown, that x(t) and y{t) satisfy the mean field 
equations presented in the previous section Plugging 
the expansion into the Master equation and substituting the 
probability distribution of the fluctuations, Yl(^^,^y,t), for 
P(n V, fty, t), we can expand the right hand side of Eq. |6] in 
powers of A^"'^^ and collect terms of the same order. From 
the leading order terms in the expansion we regain the mean- 
field equations given by Eqs. |2]and|3] The next highest order 
terms give a linear Fokker-Planck equation for distribution of 
fluctuations, which is known to have a Gaussian solution. The 
Fokker-Planck equation can be equivalently described by the 
linear Langevin equation 



(12) 



Here we have written the Langevin equation in vector form 
where ^ - (^^ , ^y), K is the drift matrix and is exactly equal to 
the Jacobian matrix from the linear stability analysis of Eqs. |2] 
and|3] and f(/) = {fx,fy) is Gaussian white noise described by 



{fi(t)fj(t')}^2Dij6(t-t') 



(13) 
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FIG. 4: Power spectrum of the stochastic Brusselator averaged 
over 500 reahzations with b = 1.9 and c = 1.0 for system sizes 
A' = lO^^(blue), A' = 10''(red), and N = lO^(black) against the power 
spectrum calculated from the van Kampen expansion of the Master 
equation(dashed) . 



The matrices K and D depend upon the mean field solutions, 
x{t) and y{t). Evaluated at the fixed point, these are by: 



K = 



D 



b-l 



l+b -b 



(14) 



(15) 



From here one can easily obtain an expression for thepower 
spectra of fluctuations of X and Y, as was shown in yfl. The 
power spectra are given by 



2{{\+bW + c^) 



(c-w2)2 + (1 +c-bfa? 



2b{{u? + \+ b) 



(C-W2)2 + (1 + C-/7)W 



(16) 



(17) 



The power spectra calculated from the van Kampen approx- 
imation are peaked at non-zero frequencies for parameters in 
the fixed point regime of the mean-field Brusselator model. 
Here the system acts analogously to a damped driven pen- 
dulum. The intrinsic fluctuations drive the system at all fre- 
quencies, thus exciting the systems natural frequency. These 
oscillations are the quasicycles originally reported in . 



m. NUMERICAL RESULTS 

We studied the intrinsically noisy Brusselator by numeri- 
cally solving the Master equation, Eq. |6] using Gillespie's 
direct method [is!]. Our goal was to compare the simulated 



FIG. 5: Power spectrum of the stochastic Brusselator averaged over 
500 realizations with b = 9.9 and c = 9.0 for system sizes A' = 
103(blue), N = 10-*(red), and N = lO'(black). 



dynamics with the van Kampen predictions as the system size 
was changed, and tuning the system through the mean-field 
bifurcation point. 

We began by looking at the case with no time scale sepa- 
ration, setting c - 1.0 as was studied in Q]. Boland showed 
that the intrinsic fluctuations give rise to quasi-cycles near the 
mean-field bifurcation point be- These quasi-cycles are char- 
acterzed by their distinct, Lorentzian-like power spectrum. In 
Fig. ID we show the power spectrum at a distance 6 - -0.1 
from the bifurcation point for three different system sizes, 
along with the van Kampen prediction. Each power spectrum 
was calculated from runs with a total time T - 5000 and av- 
eraged over 500 realizations. At large A^ the calculated power 
spectrum nicely matches the prediction, but as A^ is decreased 
the power spectrum develops a harmonic peak at 2u)q. Sim- 
ilarly, the power spectrum develops this harmonic peak as 5 
approaches 0, and the quasi-cycle becomes a noisy limit cy- 
cle iH. Note that although the van Kampen calculated power 
spectrum appears to decay more rapidly, all of the spectra do 
fall off with the same power near the peak. This can be seen 
more clearly by normalizing the peak heights. 

As c is increased, decreasing the characteristic timescale 
of y, the power spectrum undergoes a continuous change in 
shape dependent on the system size. From the van Kampen 
expansion, we expect the peak frequency to go as c'^^, how- 
ever, for smaller system sizes a peak emerges near o) - \. This 
can be seen in Fig. |5] where we plot the power spectra for the 
same A^ and 6 as were plotted in Fig. |4] Here, at N = 10^, 
the power spectrum still matches closely to the van Kampen 
approximation, while atN = 10^ it is dominated by this w = 1 
peak. 

We can begin to understand this qualitative change in be- 
havior by looking at the simulated phase space trajectories 
in the two different regimes, shown in Fig. |6] Without the 
timescale separation, the system will closely orbit the fixed 
point, shown as a red dot, where the size and shape of the 
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FIG. 6: Phase portraits for system size N = 10' at c = l.O,h = 
1.90 showing small amplitude quasicycles (top) and c = 9.0, h = 
9.90 showing both small amplitude quasicycles and large amplitude 
excursions (bottom). 



orbit can be determined from the stead-state solution to the 
Fokker-Planck equation. The covariance of these fluctuations 
can be shown to depend on both b and c, and have an overall 
scaling proportional to N^^^~. When the timescale separation 
is introduced, the same small amplitude quasi-cycles occur 
around the fixed points, but there is also now a larger, almost 
deterministic periodic trajectory around the fixed point. In this 
regime, the system will wander around the fixed point until it 
hits the edge of the linearly stable region, at which point it 
will go on a large amplitude excursion before returning again 
to the linearly stable region. The system may also undergo 
multiple large amplitude cycles before returning. Compar- 
ing the stochastic phase space trajectories with the mean field 
ones shown in Fig. |2l the difference in behavior can be traced 
back to the appearance of a line of fixed points in the c ^ oo 
limit. In that limit, longitudinal fluctuations along the degen- 
erate nullclines would overwhelm any transverse fluctuations, 
causing the system to execute the large amplitude excursions 
even deep in the fixed pint regime of the mean field phase di- 
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FIG. 7: Typical time series displaying mixed mode oscillations for 
the system with parameters N = 10^, h = 9.90, and c = 9.0. and 
Tj denote the first passage times from large to small amplitude and 
small to large amplitude oscillations respectively. 



agram. For large but finite c, transverse fluctuations result in 
the system entering the small amplitude quasicycle. In this 
regime, therefore, we expect to see mixed mode oscillations 
with switching between the small amplitude quasi cycles and 
the large amplitude excursions. 

To investigate this behavior, we consider the transitions be- 
tween oscillations of different types to be first passage pro- 
cesses. We define r , to be the first passage time going from 
small to large amplitude oscillations, and as the first pas- 
sage time from large to small. An illustration of these times is 
shown on the time series plot in Fig [7] We then collected first 
passage times at 5 different system sizes between = lO-' and 

= 10"* at a range of b values near the bifurcation point for 
c = 9.0. The first passage time distribution for Tj was found 
to be exponential, while the distribution for tl consisted of 
sharp peaks at integer multiples of the large amplitude period 
Tl that decayed within an exponential envelope. In this case 
we define the mean first passage time by binning the times 
around each peak and using (t^) - 2„ nTip„, where p„ is the 
weight of each peak. 

In Fig. [8]and Fig. |9]we plot the average first passage times 
as a function of 6 at different system sizes for t, and ti re- 
spectively. As expected, (t J decreases monotonically, while 
{ti) is monotonically increasing. The inflection points on the 
{ts) for larger systems is due to the first passage times being 
of the same order as the total simulation time. 

We then looked for a scaling function of the form (r) — 
N^^f{6N") to see how the first passage times depend on the 
system size, or noise strength. Figures|8]and|9]show the results 
of this scaling, where we found as - Q.%2,fis - -1.12 for {r j 
and ai - 0.2, /3l = 0.2 for (tl). So while both first passage 
times do appear to have a systematic dependence on the noise 
strength, it is different for each process. This is also apparent 
when (ts) and (tl) are plotted against one another as in Fig. 
[TOl As the system size increases, we see that the point where 
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FIG. 8: Average small-to-large first passage times for system sizes 
Af = (10^(A),2x 10^(n),5x 10^(o),7x lO^(V)) as a function of the 
distance from the bifurcation, S, unsealed (inset) and scaled by the 
system size with exponents a = 0.82 and yS = 1.12. 



N= 10000c = 9 




FIG. 9: Average large-to-small first passage times for system sizes 
7^ = {103(A), 2 X 10^(0), 5 X 10^(0), 7 X 10''(V)) as a function of the 
distance from the bifurcation, <5, unsealed (inset) and scaled by the 
system size with exponents a = = 0.20. 




FIG. 10: Average first passage times as a function of S for 5 difi^erent 
systems sizes(A' = 10\2 x 10\5 x 10^7 x 10^ lO"*). As the sys- 
tem size increases the point at which (t,) = (tl) shifts to the left, 
approaching the point where the amplitude jumps in the mean-field 
singular Hopf bifurcation. 



small amplitude oscillations dominating far away from the bi- 
furcation and the large amplitude oscillations dominating for 
6 > 6q, where the mean field amplitude shows a jump due 
to the singular nature of the Hopf bifurcation ifTTl . For finite 
c, the scaling shown in Fig. [8] is expected to breakdown at a 
large enough value of A^, where the Van Kampen expansion 
becomes a good approximation. This breakdown occurs at 
larger and larger A^s as c is increased and based on the nature 
of the mean field phase portraits, we expect that as c ^ oo, 
the mixed mode oscillations will persist for all system sizes. 



the two timescales cross shifts to the point where we see an 
amplitude jump in the mean-field system shown in Fig. [3] 
Beyond this value of 6 the system becomes dominated by the 
large amplitude oscillations. It is also interesting to note that 
this crossing point occurs at the same (r) for all of the system 
sizes we analyzed. The scaling relations for (tl), and (tj), 
imply that the scaling functions f{x) decay as x^^" in order for 
the system to have a well-defined large limit. In this limit, 
(tl) X S, and (tJ r; These predictions are consistent 

with the numerical results in Fig. |8] 

In the simulations for large c, there is no signature of the 
Hopf bifurcation in the noisy system. Instead, there is switch- 
ing between small and large amplitude oscillations, with the 



IV. DISCUSSION 

In this work we analyzed the eff'ects of demographic or in- 
trinsic noise in the zero dimensional Brusselator, a prototypi- 
cal chemical oscillator representing a well-mixed system, and 
showed that noise can qualitatively change the nature of the 
system in the limit of large timescale separation between the 
inhibitor and the activator We verified previous work showing 
that the noise induced quasi-cycles become limit cycle oscil- 
lations as the system is tuned past the bifurcations point and 
then showed that this transition occurs earlier for smaller sys- 
tems. We then looked at what happens when a time scale sep- 
aration between the two variable is introduced and found that 
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a noise dependent switching between small and large ampli- 
tude oscillation occurs. Based on the scaling relations for the 
first passage times, we speculate that carefully taking the lim- 
its —» oo and c — > oo simultaneously will lead to a new set 
of deterministic equations. This will be the subject of future 
work. 

Previous analysis of the Brusselator in the limit where the 
activator time scale is much shorter than the inhibitor time 
scale has shown the presence of stochastic resonance lfl2l . 
Noise-induced mixed mode oscillations in systems whose 
mean field equations exhibit a singular Hopf bifurcation has 
also been demonstrated in earlier work H. Such studies 
demonstrate the singular nature of noise close to the bifur- 
cation. Our analysis of the Brusselator in the regime where 
the inhibitor time scale is fast indicates a singular effect of 
intrinsic noise in regions far away from the bifurcation. The 
origin of noise-induced mixed mode oscillations in regimes 
where the fixed point should be stable is the emergence of a 
line of fixed points in the c — > oo limit. For such a line of 
fixed points, noise has a large destabilizing influence simi- 



lar to "soft modes" in classical statistical mechanics systems 
where fluctuations can destroy a fixed point |[l9l . 

We are currently working on extending the work presented 
here by looking at how noise-induced bistability affects the 
dynamics of the spatially extended Brusselator system, where 
diffusion is relevant. Our approach is to look at a lattice of 
diffusively couple point oscillator and apply the same tech- 
niques presented here. Preliminary simulations have shown 
similar switching behaviors, but with different first passage 
time distributions. 
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